# Get data from NOAA database
url = "ftp://ftp.cdc.noaa.gov/Datasets/COBE2/sst.mon.mean.nc"
bin = getBinaryURL(url)
writeBin(bin, "mon.sst.nc")
SSTData <- nc_open(filename = "mon.sst.nc")
rm(bin)
file.remove("mon.sst.nc")

# Set variables from dataset
lon <- ncvar_get(SSTData, varid = "lon")
lat <- ncvar_get(SSTData, varid = "lat")
time <- ncvar_get(SSTData, varid = "time")
sst <- ncvar_get(SSTData, varid = "sst")

# Set time as a date to match extracted data
dim(sst)
SSTData$dim$time$units
sst.date <- as.Date("1891-01-01") + time

# First define a helper function that extracts middle value (or average of 2 middle values) from a vector).
midval <- function(x) {
  if(length(x)%%2 == 1) {
    return(x[ceiling(length(x)/2)])
  } else return(mean(c(x[length(x)/2], x[length(x)/2 + 1]), na.rm = T))
  
}

# Define a function for extracting mean SST value at a particular location from a particular time period
sst_extract <- function(start_date, end_date, data_source,
                        lon_index, lat_index,
                        summer = F,
                        fun = mean, ...) {
  require("lubridate")
  if(summer) {
    if(lat_index <= 90) {
      month(start_date) <- 6
      end_date <- start_date
      month(end_date) <- 8
    } else if(lat_index > 90) {
      month(start_date) <- 12
      end_date <- start_date
      month(end_date) <- 2
      year(end_date) <- year(end_date) + 1
    }
  }
  
  start_date_units <- as.numeric(start_date - as.Date("1891-01-01"))
  end_date_units <- as.numeric(end_date - as.Date("1891-01-01"))
  start_date_index <- which(time == start_date_units)
  end_date_index <- which(time == end_date_units)
  
  # ...and extract data
  sst_extract <- data_source[lon_index, lat_index, seq(start_date_index, end_date_index)]
  # Here we assign data for i-th record to our list...
  
  return(fun(sst_extract, ...))
}


# Create function which finds the closest coordinate in SST dataset from our extracted dataset
find_close <- function(datasource, lon_ix, lat_ix, step = 3, summer = F, start_date, end_date, diagn = F) {
  lon_steps = seq(lon_ix - step, lon_ix + step)
  lat_steps = seq(lat_ix - step, lat_ix + step)
  
  longitudes = numeric(length(lon_steps)) + 1
  latitudes = numeric(length(lat_steps)) + 1
  names(longitudes) = lon_steps
  names(latitudes) = lat_steps
  
  for(i in 1:length(longitudes)) {
    
    sst_vals = sst_extract(start_date = start_date, end_date = end_date, data_source = datasource,
                           lat_index = lat_ix,
                           lon_index = as.numeric(names(longitudes)[i]),
                           summer = summer,
                           fun = function(x) return(x))
    longitudes[i] = any(!is.na(sst_vals))
    
  }
  
  for(i in 1:length(latitudes)) {
    
    sst_vals = sst_extract(start_date = start_date, end_date = end_date, data_source = datasource,
                           lon_index = lon_ix,
                           lat_index = as.numeric(names(latitudes)[i]),
                           summer = summer,
                           fun = function(x) return(x))
    latitudes[i] = any(!is.na(sst_vals))
    
  }
  
  if (diagn) {
    return(rbind(longitudes, latitudes))
  } else {
    
    if(any(longitudes == 1)) {
      for (k in 1:step) {
        if (longitudes[step + 1 + k] == 1) {
          lon_new = as.numeric(names(longitudes)[step + 1 + k])
          break
        } else if (longitudes[step + 1 - k] == 1) {
          lon_new = as.numeric(names(longitudes)[step + 1 - k])
          break
        }
      }
      lon_ix = lon_new
    } else {
      
      for (k in 1:step) {
        if (latitudes[step + 1 + k] == 1) {
          lat_new = as.numeric(names(latitudes)[step + 1 + k])
          break
        } else if (latitudes[step + 1 - k] == 1) {
          lat_new = as.numeric(names(latitudes)[step + 1 - k])
          break
        }
      }
      lat_ix = lat_new
      
    }
    
    
    return(c(lon_ix, lat_ix))
    
  }
}

# Create a function to identify survey season from sampling period
season_extract <- function(start_date, end_date, lat_index) {
  require('lubridate')
  
  N_seasons <- c(rep('win', 2), rep('spr', 3), rep('sum', 3), rep('aut', 3), 'win')
  S_seasons <- c(rep('sum', 2), rep('aut', 3), rep('win', 3), rep('spr', 3), 'sum')
  
  if (as.numeric(difftime(end_date, start_date)) < 100) {
    
    if (lat_index <= 90) {
      return(names(sort(table(N_seasons[month(start_date):month(end_date)]), decreasing = T))[1])
    } else if (lat_index > 90) {
      return(names(sort(table(S_seasons[month(start_date):month(end_date)]), decreasing = T))[1])
    }
  } else return('multi')
}


# Assign month names to number values
months <- c(Jan = "01", Feb = "02", Mar = "03",
            Apr = "04", May = "05", Jun = "06",
            Jul = "07", Aug = "08", Sep = "09",
            Oct = "10", Nov = "11", Dec = "12")

# Set extracted data as dataframe
percentonly <- as.data.frame(percentonly)

# Create columns for values to go into in dataframe
percentonly$average_SST <- NA
percentonly$middle_SST <- NA
percentonly$sd_SST <- NA
percentonly$start_month <- NA
percentonly$end_month <- NA
percentonly$average_SST_summer <- NA
percentonly$middle_SST_summer <- NA
percentonly$sd_SST_summer <- NA

# Allow coordinates to round to the nearest 0.5 degree to best match to SST dataset
coord_grid <- function(coord) floor(coord) + 0.5

# Rename in case any variables accidentally get changed, so original is saved
percentonly_ <- percentonly

for (i in 1:nrow(percentonly_)) {
  # cat(i); cat("\n")
  if (grepl("^[A-Z]{1}[a-z]{2}$",
            percentonly_[i, "Month"])) {
    # this condition looks for cases with one month
    
    start.month <- end.month <- percentonly_[i, "Month"]
    start.yr <- end.yr <- percentonly_[i, "Year"]
    
    
  } else if (grepl("^([A-Z]{1}[a-z]{2}), (?:[A-Z]{1}[a-z]{2}, )*([A-Z]{1}[a-z]{2})$",
                   percentonly_[i, "Month"])) {
    # this condition looks for lists of month separated by ", "
    
    start.month <- gsub("^([A-Z]{1}[a-z]{2}), (?:[A-Z]{1}[a-z]{2}, )*([A-Z]{1}[a-z]{2})$",
                        replacement = "\\1",
                        percentonly_[i, "Month"])
    end.month <- gsub("^([A-Z]{1}[a-z]{2}), (?:[A-Z]{1}[a-z]{2}, )*([A-Z]{1}[a-z]{2})$",
                      replacement = "\\2",
                      percentonly_[i, "Month"])
    start.yr <- percentonly_[i, "Year"]
    end.yr <- ifelse(which(names(months) == end.month) > which(names(months) == start.month),
                     start.yr, start.yr + 1)
    
  } else if (grepl("^([A-Z]{1}[a-z]{2})-[A-Z]{1}[a-z]{2}, (?:[A-Z]{1}[a-z]{2}-[A-Z]{1}[a-z]{2}, )*[A-Z]{1}[a-z]{2}-([A-Z]{1}[a-z]{2})$",
                   percentonly_[i, "Month"])) {
    # this condition looks for ranges of months signified by "-" in multiple ranges
    
    start.month <- gsub("^([A-Z]{1}[a-z]{2})-[A-Z]{1}[a-z]{2}, (?:[A-Z]{1}[a-z]{2}-[A-Z]{1}[a-z]{2}, )*[A-Z]{1}[a-z]{2}-([A-Z]{1}[a-z]{2})$",
                        replacement = "\\1",
                        percentonly_[i, "Month"])
    end.month <- gsub("^([A-Z]{1}[a-z]{2})-[A-Z]{1}[a-z]{2}, (?:[A-Z]{1}[a-z]{2}-[A-Z]{1}[a-z]{2}, )*[A-Z]{1}[a-z]{2}-([A-Z]{1}[a-z]{2})$",
                      replacement = "\\2",
                      percentonly_[i, "Month"])
    start.yr <- percentonly_[i, "Year"]
    end.yr <- ifelse(which(names(months) == end.month) > which(names(months) == start.month),
                     start.yr, start.yr + 1)
    
  } else if (grepl("^([A-Z]{1}[a-z]{2})-([A-Z]{1}[a-z]{2})$",
                   percentonly_[i, "Month"])) {
    # this condition looks for ranges of months signified by "-"
    
    start.month <- gsub("^([A-Z]{1}[a-z]{2})-([A-Z]{1}[a-z]{2})$",
                        replacement = "\\1",
                        percentonly_[i, "Month"])
    end.month <- gsub("^([A-Z]{1}[a-z]{2})-([A-Z]{1}[a-z]{2})$",
                      replacement = "\\2",
                      percentonly_[i, "Month"])
    start.yr <- percentonly_[i, "Year"]
    end.yr <- ifelse(which(names(months) == end.month) > which(names(months) == start.month),
                     start.yr, start.yr + 1)
    
  } else if (percentonly_[i, "Month"] == 0) {
    start.yr <- end.yr <- percentonly_[i, "Year"]
    if(percentonly_[i, "Lat"] > 0) start.month <- end.month <- "Jul"
    if(percentonly_[i, "Lat"] < 0) start.month <- end.month <- "Jan"
    
  } else start.month <- end.month <- middle.month <- -999
  
  if (start.month!= -999) {
    lat_index <- which(lat == coord_grid(percentonly_[i, "Lat"])) # ...extract relevant indexes...
    lon_index <- which(lon == coord_grid(ifelse(percentonly_[i, "Lon"] < 0,
                                                360 + percentonly_[i, "Lon"],
                                                percentonly_[i, "Lon"])))
    
    start_date <- as.Date(paste(start.yr, "-", months[start.month], "-01", sep = ""))
    end_date <- as.Date(paste(end.yr, "-", months[end.month], "-01", sep = ""))
    
    
    if(any(is.na(sst_extract(start_date, end_date, sst,
                             lon_index, lat_index, summer = F, function(x) return (x))))) {
      lon_index = find_close(datasource = sst, lon_ix = lon_index, lat_ix = lat_index,
                             step = 10, summer = F, start_date = start_date, end_date = end_date)[1]
      lat_index = find_close(datasource = sst, lon_ix = lon_index, lat_ix = lat_index,
                             step = 10, summer = F, start_date = start_date, end_date = end_date)[2]
    }
    
    percentonly_[i, "start_month"] <- as.numeric(format(start_date, "%m"))
    percentonly_[i, "end_month"] <- as.numeric(format(end_date, "%m"))
    
    percentonly_[i, "average_SST"] <- sst_extract(start_date, end_date, sst,
                                                  lon_index, lat_index, summer = F, mean, na.rm = T)
    percentonly_[i, "middle_SST"] <- sst_extract(start_date, end_date, sst,
                                                 lon_index, lat_index, summer = F, midval)
    percentonly_[i, "sd_SST"] <- sst_extract(start_date, end_date, sst,
                                             lon_index, lat_index, summer = F, sd, na.rm = T)
    percentonly_[i, "average_SST_summer"] <- sst_extract(start_date, end_date, sst,
                                                         lon_index, lat_index, summer = T, mean, na.rm = T)
    percentonly_[i, "middle_SST_summer"] <- sst_extract(start_date, end_date, sst,
                                                        lon_index, lat_index, summer = T, midval)
    percentonly_[i, "sd_SST_summer"] <- sst_extract(start_date, end_date, sst,
                                                    lon_index, lat_index, summer = T, sd, na.rm = T)
    percentonly_[i, "season"] <- season_extract(start_date, end_date, lat_index)
    
    
  } else if (start.month == -999) {
    
  } else stop("Error in month formatting\n")
}



# rename back to original once finished
percentonly <- percentonly_

## Check sst data coverage
out <- array(0, dim = c(360, 180, 2040))
plot(1:360, 1:360, type = "n", ylim = c(1,180))
for (i in 1:360) {
  for (j in 1:180) {
    for(k in 1:2040) {
      if(!is.na(sst[i,j,k])) {
        # points(lon[i], lat[j])
        out[i,j,k] <- out[i,j,k] + 1
      }
    }
  }
}

out1 <- rowSums(out, dims = 2)
dim(out1)
plot(1:360, 1:360, type = "n", ylim = c(180,1))
for (i in 1:360) {
  for(j in 1:180) {
    points(i, j, cex = 0.001 * out1[i,j])
  }
}

# Check if averages make sense with middle month values
ggplot(percentonly, aes(x = average_SST, y = middle_SST)) + geom_point()

# remove NaN results
TStemp <- subset(percentonly, (middle_SST == "NaN"))

# Check correlation between year and average SST
cor.test(percentonly$Year, percentonly$average_SST)
ggplot(percentonly, mapping = aes(Year, average_SST, col = abs(Lat))) + 
  geom_jitter() +
  geom_smooth()


# Check by lat values without <25
ggplot(percentonly, mapping = aes(Year, middle_SST, col = abs(Lat))) + 
  geom_jitter() +
  geom_smooth()
ggplot(percentonly, mapping = aes(Year, average_SST, col = abs(Lat))) +
  geom_point() +
  geom_smooth(method = "lm")
ggplot(percentonly, mapping = aes(Year, average_SST_summer, col = abs(Lat))) +
  geom_point() +
  geom_smooth(method = "lm")